#ifndef _VECTORFUNCTION_NEW_H_
#define _VECTORFUNCTION_NEW_H_

#include <iostream>
#include <cmath>
#include <ctime>
#include "REAL.H"
#include "Box.H"
#include "Vector.H"
#include "RealVect.H"
#include "FArrayBox.H"
#include "BoxIterator.H"
#include "BoxLayoutData.H"
#include "DisjointBoxLayout.H"
#include "LevelData.H"
#include "DataIterator.H"
#include "IntVectSet.H"
#include "LayoutIterator.H"
#include "UsingNamespace.H"

#include "NamespaceHeader.H"

enum FABType
{
    ModeCellAverage = 1,
    ModeCellCenter = 2
};

class VectorFunction
{
protected:
    FABType FABMode = ModeCellAverage;

public:
    VectorFunction(){};
    virtual ~VectorFunction(){};
    /// compute function value for a single point (x,t)
    virtual const Vector<Real> operator()(const RealVect &x, const Real &time = 0) const = 0;
    /// Return the number of components of the returned vector.
    virtual const int getNumComponents(void) const { return 0; };
    /// Set the values of a FArrayBox.
    virtual void setFAB(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time = 0, const int startComp = 0) const;
    /// Set LevelData to averaged values of the vector function on its boxes.
    virtual void setLFA(LevelData<FArrayBox> &lfa, const RealVect &origin, const Real &dx, const Real &time = 0, const int startComp = 0) const;
    /// ModeCellAverge = 1, ModeCellCenter = 2
    void setFABmode(const int FABModeCode) { FABMode = (enum FABType)(FABModeCode); };
    virtual void setFABCellCenter(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time = 0, const int startComp = 0) const;
    virtual void setFABCellAverage(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time = 0, const int startComp = 0) const;
    /// Set the Laplacian values of a FArraybox (we combine Origins and Lapacian-results by storing Lapacian after Origins at indexStore).
    virtual void setLaplacian(FArrayBox &fab, const RealVect &origin, const Real &dx, const int startComp, const int indexStore) const;
};

void VectorFunction::setFAB(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time, const int startComp) const
{
    switch ((this->FABMode))
    {
    case ModeCellAverage:
        setFABCellAverage(fab, origin, dx, time, startComp);
        break;
    case ModeCellCenter:
        setFABCellCenter(fab, origin, dx, time, startComp);
        break;
    default:
        break;
    }
}

void VectorFunction::setFABCellCenter(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time, const int startComp) const
{
    Box b = fab.box();
    BoxIterator bit(b);
    const int NumComponents = this->getNumComponents();
    for (bit.begin(); bit.ok(); ++bit)
    {
        IntVect iv = bit();
        RealVect temp = iv;
        Vector<Real> valCellCenter = (*this)((temp + 0.5) * dx + origin, time);
        for (int i = 0; i < NumComponents; i++)
            fab.set(iv, startComp + i, valCellCenter[i]);
    }
}

void VectorFunction::setFABCellAverage(FArrayBox &fab, const RealVect &origin, const Real &dx, const Real &time, const int startComp) const
{
    Box b = fab.box();
    BoxIterator Bit(b);
    const int NumComponents = this->getNumComponents();
    const int InteralNum = 3;
    Vector<Real> SimpsonRatio;
    SimpsonRatio.push_back(1.0 / 6);
    SimpsonRatio.push_back(2.0 / 3);
    SimpsonRatio.push_back(1.0 / 6);
/// an ugly way to implement numerical integral both on DIM= 2, 3.
#if CH_SPACEDIM == 2
    for (Bit.begin(); Bit.ok(); ++Bit)
    {
        IntVect iv = Bit();
        Vector<Real> rvx;
        Vector<Real> rvy;
        Vector<Real> valCellAverage;
        for (int i = 0; i < InteralNum; i++)
        {
            rvx.push_back((iv[0] + i * 1.0 / (InteralNum - 1)) * dx + origin[0]);
            rvy.push_back((iv[1] + i * 1.0 / (InteralNum - 1)) * dx + origin[1]);
        }
        for (int k = 0; k < NumComponents; k++)
            valCellAverage.push_back(0);
        for (int i = 0; i < InteralNum; i++)
        {
            for (int j = 0; j < InteralNum; j++)
            {
                RealVect rv(rvx[i], rvy[j]);
                Vector<Real> temp = (*this)(rv, time);
                for (int k = 0; k < NumComponents; k++)
                    valCellAverage[k] += SimpsonRatio[i] * SimpsonRatio[j] * temp[k];
            }
        }
        for (int i = 0; i < NumComponents; i++)
            fab.set(iv, startComp + i, valCellAverage[i]);
    }
#elif CH_SPACEDIM == 3
    for (Bit.begin(); Bit.ok(); ++Bit)
    {
        IntVect iv = Bit();
        Vector<Real> rvx;
        Vector<Real> rvy;
        Vector<Real> rvz;
        Vector<Real> valCellAverage;
        for (int i = 0; i < InteralNum; i++)
        {
            rvx.push_back((iv[0] + i * 1.0 / (InteralNum - 1)) * dx + origin[0]);
            rvy.push_back((iv[1] + i * 1.0 / (InteralNum - 1)) * dx + origin[1]);
            rvz.push_back((iv[2] + i * 1.0 / (InteralNum - 1)) * dx + origin[2]);
        }
        for (int k = 0; k < NumComponents; k++)
            valCellAverage.push_back(0);
        for (int i = 0; i < InteralNum; i++)
        {
            for (int j = 0; j < InteralNum; j++)
            {

                for (int l = 0; l < InteralNum; l++)
                {
                    RealVect rv(rvx[i], rvy[j], rvz[l]);
                    Vector<Real> temp = (*this)(rv, time);
                    for (int k = 0; k < NumComponents; k++)
                        valCellAverage[k] += SimpsonRatio[i] * SimpsonRatio[j] * SimpsonRatio[l] * temp[k];
                }
            }
        }
        for (int i = 0; i < NumComponents; i++)
            fab.set(iv, startComp + i, valCellAverage[i]);
    }
#else
    DIM > 3 undefined !
#endif
}

void VectorFunction::setLFA(LevelData<FArrayBox> &lfa, const RealVect &origin, const Real &dx, const Real &time, const int startComp) const
{
    DataIterator dir = lfa.dataIterator();
    for (dir.begin(); dir.ok(); ++dir)
        this->setFAB(lfa[dir], origin, dx, time, startComp);
}

void VectorFunction::setLaplacian(FArrayBox &fab, const RealVect &origin, const Real &dx, const int startComp, const int indexStore) const
{
    const int NumComponents = this->getNumComponents();
    Box b = fab.box();
    const int *bloend = b.loVect();
    const int *bhiend = b.hiVect();
    BoxIterator bit(b);
    for (bit.begin(); bit.ok(); ++bit)
    {
        IntVect iv = bit();
        Real temp;
        int flag = 1;
        for (int j = 0; j < CH_SPACEDIM; j++)
            flag *= ((iv[j] - bloend[j]) * (iv[j] - bloend[j] - 1) * (iv[j] - bhiend[j]) * (iv[j] - bhiend[j] + 1));
        if (flag)
        {
            for (int i = 0; i < NumComponents; i++)
            {
                temp = 0;
                for (int j = 0; j < CH_SPACEDIM; j++)
                {
                    temp += (-fab(iv + 2 * BASISV(j), startComp + i) + 16 * fab(iv + BASISV(j), startComp + i) - 30 * fab(iv, startComp + i) + 16 * fab(iv - BASISV(j), startComp + i) - fab(iv - 2 * BASISV(j), startComp + i));
                }
                temp /= (12 * dx * dx);
                fab.set(iv, indexStore + i, temp);
            }
        }
        else
        {
            for (int i = 0; i < NumComponents; i++)
            {
                fab.set(iv, indexStore + i, 0);
            }
        };
    }

}

#include "NamespaceFooter.H"

#else
// do nothing
#endif